Cite this Article
Lv Xinjing, Xu Xinyu, Feng Qi, Zhang Bin, Ding Yingchun, Liu Qiang. High-contrast imaging based on wavefront shaping to improve low signal-to-noise ratio photoacoustic signals using superpixel method. Chinese Physics B, 2020, 29(3): 034301  
Permissions
High-contrast imaging based on wavefront shaping to improve low signal-to-noise ratio photoacoustic signals using superpixel method
Lv Xinjing1, Xu Xinyu1, Feng Qi1, Zhang Bin2, Ding Yingchun1, †, Liu Qiang2, ‡
College of Mathematics and Physics, Beijing University of Chemical Technology, Beijing 100029, China
State Key Laboratory of Precision Measurement Technology and Instruments, Department of Precision Instruments, Tsinghua University, Beijing 100084, China

 

† Corresponding author. E-mail: dingyc@mail.buct.edu.cn qiangliu@mail.tsinghua.edu.cn

Project supported by the National Key Research and Development Program of China (Grant No. 2017YFB1104500), the Beijing Natural Science Foundation, China (Grant No. 7182091), the National Natural Science Foundation of China (Grant No. 21627813), and the Research Projects on Biomedical Transformation of China-Japan Friendship Hospital (Grant No. PYBZ1801).

Abstract

Photoacoustic (PA) imaging has drawn tremendous research interest for various applications in biomedicine and experienced exponential growth over the past decade. Since the scattering effect of biological tissue on ultrasound is two- to three-orders magnitude weaker than that of light, photoacoustic imaging can effectively improve the imaging depth. However, as the depth of imaging further increases, the incident light is seriously affected by scattering that the generated photoacoustic signal is very weak and the signal-to-noise ratio (SNR) is quite low. Low SNR signals can reduce imaging quality and even cause imaging failure. In this paper, we proposed a new wavefront shaping and imaging method of low SNR photoacoustic signal using digital micromirror device (DMD) based superpixel method. We combined the superpixel method with DMD to modulate the phase and amplitude of the incident light, and the genetic algorithm (GA) was used as the wavefront shaping algorithm. The enhancement of the photoacoustic signal reached 10.46. Then we performed scanning imaging by moving the absorber with the translation stage. A clear image with contrast of 8.57 was obtained while imaging with original photoacoustic signals could not be achieved. The proposed method opens new perspectives for imaging with weak photoacoustic signals.

PACS: ;43.35.Ud;;42.25.Dd;;42.25.Fx;;42.30.-;
Keyword:;photoacoustic imaging;;wavefront shaping;;superpixel method;;high contrast imaging;
1. Introduction

Noninvasive biomedical imaging plays an important role in modern medicine. Imaging of blood vessel, brain tissue, and other biological tissues is of great help to disease diagnosis, disease treatment, and drug delivery.[1–3] When light propagates in biological tissues, it undergoes multiple scattering and rapidly becomes diffused light.[4] Usually, random speckles are formed due to multiple scattering when light passes through the scattering media. Traditional optical imaging methods such as confocal microscopy,[5] optical coherence tomography (OCT),[6] and two- or three-photon microscopy[7] utilize ballistic photons to image and often have high resolution. But influenced by scattering, the imaging depth of these methods is limited, usually less than one transport mean free path.[8] Thus, they can only be applied to superficial biological tissues.

The scattering effect of biological tissues on ultrasound is two- to three-orders magnitude weaker than that of light.[9] Based on photoacoustic effect, photoacoustic imaging[10,11] combines the advantages of light absorption contrast and ultrasound transmission. Photoacoustic effect refers to the phenomenon that the absorber absorbs the energy of the incident laser and excites ultrasound signal. When the laser irradiates the absorber, the absorber absorbs the light energy of the incident light and converts it into heat, causing local temperature rise and thermal expansion. As a result, a raised pressure is generated in the area of the absorption body irradiated by the light. If the laser pulse width is shorter than the stress relaxation time and the thermal relaxation time, heat and pressure will not diffuse to the surrounding area apparently within the duration of the laser pulse. With thermal diffusion, the absorber produces thermal-elastic deformation and transmits pressure outward which forms ultrasound. The generated ultrasound signal is the photoacoustic signal. The amplitude of the photoacoustic signal is determined by the laser fluence F and the absorption coefficient μa of the absorber

where Γ denotes the Grüneisen coefficient of the absorber. Optical absorption contrast image of the absorber can be reconstructed by collecting the photoacoustic signals. In addition to high resolution and contrast, photoacoustic imaging has the advantage of imaging depth since ultrasound is less affected by scattering.

Although ultrasound can travel farther than light in the scattering media, the incident laser beam is still seriously affected by scattering. As the imaging depth increases, the influence of scattering on the incident light is aggravated and the generated photoacoustic signal weakens or even disappears. To enhance the amplitude of the photoacoustic signal, it is necessary to reduce the effect of scattering on the incident light and focus light as much as possible to the target point.

Wavefront shaping methods[12] have been developed to focus light through scattering media noninvasively. The basic principle is to modulate the wavefront phase and amplitude of the incident wave through an optical modulator to compensate for the light distortion caused by scattering. Wavefront shaping methods including transmission matrix method,[13–15] iterative optimization method,[16] optical phase conjugation (OPC),[17,18] and time-reversed ultrasonically encoded (TRUE) optical focusing technique[19] have made remarkable progress, especially in resolution and imaging speed. The light intensity of the target point is detected as the feedback signal for iterative optimization wavefront shaping. Focusing light through the scattering media by wavefront shaping has drawn much attention since the light detector cannot detect the light intensity of the target point which is inside the scattering media noninvasively. Although the TRUE method can detect the ultrasonically encoded light coming from the focal point using a light detector placed outside the scattering media, this method is hindered by complex experimental optical setup and cumbersome experimental procedures.[19]

The amplitude of the photoacoustic signal is proportional to the light fluence of the target area that the photoacoustic signal can be used as the feedback signal for wavefront shaping of the incident light.[20] Combining photoacoustic imaging and wavefront shaping further reduces the influence of scattering on the incident light and the photoacoustic signals collected for image reconstruction are significantly enhanced. In 2011, Kong et al. first proposed the idea of photoacoustic guided wavefront shaping.[21] They modulated the wavefront of the incident light with a deformable mirror (DM) to maximize the photoacoustic signal and achieved light focusing. In 2013, Chaigne et al. noninvasively measured the photoacoustic transmission matrix of the scattering media with the phase-only spatial light modulator (SLM) exploiting the photoacoustic effect.[22] Light focusing through the scattering media could be achieved utilizing the measured photoacoustic transmission matrix and the enhancement of the photoacoustic signal reached about 6. In 2013, Caravaca-Aguirre et al. combined the SLM with the genetic algorithm to iteratively optimize the wavefront of the incident light with the photoacoustic feedback, and the enhancement reached 10.[23] The SLM is commonly used in photoacoustic wavefront shaping researches to modulate the wavefront of the incident light.[24–26] The low refresh rate of the SLM, about 100 Hz, restricts its application in vivo imaging since the speckle decorrelation time[27] of biological tissue is short. Much faster refresh rate can be achieved with digital micromirror device (DMD), up to 32 kHz.[28] However, the DMD surface can only endure pulse laser fluence no more than 200 μJ/cm2 while the SLM surface can endure 40 mJ/cm2. Lower pulsed laser fluence leads to weaker photoacoustic signal and lower SNR. In 2014, Wang et al. first used DMD for photoacoustic wavefront shaping.[29] They raised the SNR to 3.9 by averaging 64 acquisitions of photoacoustic signals that the optimization process took about 2 hours with 1-kHz laser pulse repetition rate. The data averaging method is time consuming and can only slightly improve the SNR. In 2016, Tzang et al. proposed a lock-in detection method for photoacoustic signals and improved the SNR by at least one order of magnitude.[30] With 20-kHz laser pulse repetition rate, the linear photoacoustic wavefront shaping process took 6 minutes and the nonlinear photoacoustic wavefront shaping process took about 1 hour. The enhancement of the photoacoustic signals reached 9 and 16, respectively. However, this method requires complicated circuits for analog signal processing and expensive light modulation devices. In 2019, Sun et al. combined wavelet denoising with correlation detection and improved the SNR of the photoacoustic signals by 6.5 times.[31] They used DMD for photoacoustic wavefront shaping and the enhancement of the photoacoustic signal reached 7.83. The optimization process was completed in 30 minutes with only 10-Hz laser pulse repetition rate.

The effect of phase and amplitude co-modulation is much better than pure amplitude modulation.[32] As a binary amplitude modulation digital micromirror device, DMD cannot modulate the phase of the incident light directly. The superpixel method[33] has been developed to achieve amplitude and phase co-modulation with DMD. In 2018, Jia et al. first achieved light focusing through the scattering media with DMD-based superpixel method.[34] So far, the DMD-based superpixel method has never been used in photoacoustic.

In this paper, we applied the DMD-based superpixel method to photoacoustic wavefront shaping and imaging research of low SNR photoacoustic signals for the first time. We extracted the photoacoustic signals submerged in noise with wavelet denoising and correlation detection method. Then we combined the superpixel method with DMD to modulate the phase and amplitude of the incident light and used the genetic algorithm for iterative optimization of the wavefront. Enhancement of the photoacoustic signal reached 10.46. Scanning imaging was performed by moving the sample on a two-dimensional translation stage. After optimization, a clear image with contrast of 8.57 was obtained while imaging with original photoacoustic signals directly could not be achieved.

2. Experiment setup

The experiment setup is shown in Fig. 1. A 532-nm pulsed laser (LABest, SGR-10) with 10-Hz pulse repetition frequency and ∼ 800-μJ single pulse energy was used. The beam splitter divided the laser beam into two beams. The small fraction reflected light was directed to the photodiode (Hamamatsu Photonics, S5971) to monitor the pulse energy. The photoacoustic signal would be divided by the light pulse energy to eliminate the influence of the light pulse energy fluctuation, and then it can be used as the feedback signal of wavefront optimization. Most of the light transmitted through the BS and illuminated the DMD (Texas Instruments, DLP6500) after expansion and collimation. The light reflected by the DMD passed through a 4f system consisted of two lenses, one with a focal length of 150 mm and the other with a focal length of 50 mm, and a spatial filter whose place was carefully selected that only the +1 order diffraction beam could pass it through. The outgoing beam was focused on the scattering diffuser (Edmund, Diffuser 83419) through a 10× objective (NA = 0.25) before hitting the absorption sample. The absorption sample in the water tank was placed 1.6 cm behind the scattering diffuser to produce speckles, in which the size of the speckle grains was approximate to that of the absorption sample. We used black nylon wires with a diameter of 150 μm as absorption sample. A focused water-immersed ultrasonic transducer (SIUI, 5Z10SJ30DJ) with 5-MHz central frequency was used to detect the PA signals generated from the absorption sample. The frequency of the photoacoustic signal generated by the absorber is related to its size. According to the theoretical formula:[35] fUT ≈ 0.66cs/Da, where cs is the velocity of ultrasound in water and Da is the diameter of the absorber, the central frequency of the ultrasound signal generated by the black nylon wires with a diameter of 150 μm is about 6.6 MHz, which is well within the detection bandwidth of the 5-MHz ultrasonic transducer. The detected signal was amplified by an amplifier (Mini-Circuits, ZFL-500LN+) and then collected by the oscilloscope (LeCroy, 806Zi-A). The water tank was fixed on the two-dimensional translation stage (Daheng, GCD-203050M/301101M) and moved in the x–z plane to scan and image. The computer processed the signals collected by the oscilloscope in real time and refreshed the mask on the DMD according to the feedback signals.

Fig. 1. Schematic diagram of the experimental setup. BS: beam splitter; PD: photodiode; DMD: digital micromirror device; L1, L2: lens; SF: spatial filter; OBJ: objective; S: scattering diffuser; UT: ultrasonic transducer; Amp: amplifier; OSC: oscilloscope.
3. Experimental results and discussion

Figure 2(a) shows a representative random photoacoustic signal collected when a random superpixel mask was loaded on the DMD. The blue curve represents the unprocessed raw signal. Affected by white noise, clutters caused by static electricity in equipment and the noise of amplifier, the SNR of the raw PA signals was quite low. Here, the SNR is defined as the ratio of the photoacoustic signal amplitude to the standard deviation of the noise without laser irradiation.[29] We calculated the SNR of 100 random raw signals corresponding to 100 random patterns on the DMD to be 2.14. Although photoacoustic imaging is a practical method, it is difficult to directly achieve high-quality imaging with such low signal-to-noise ratio signals.

Fig. 2. (a) The random signal before and after denoising. The blue curve represents the raw signal without denoising. The red curve represents the signal after denoising. (b) Scanning image reconstructed with the raw signals. (c) Scanning image reconstructed with the denoised signals. The insets represented the random superpixel mask.

Scanning imaging was performed by moving the absorption sample in x–z plane on the two-dimensional translation stage. We scanned 60 points along the x direction and 65 points along the z direction with a scanning step of 150 μm. At each point, 10 signals were collected for data averaging. The maximum amplitude projection image of the crossed nylon wires could be reconstructed by calculating the photoacoustic signal amplitude of each point. As shown in Fig. 2(b), the absorption sample is completely unrecognizable in the image reconstructed with the unprocessed raw signals since the useful PA signals are submerged by the white noise and clutters.

In order to improve the SNR of the PA signals, we utilized wavelet analysis for signal denoising and correlation detection to remove the invalid signals in which the clutters and the photoacoustic signal overlapped.[31] Wavelet analysis[36,37] is an ideal tool for signal time –frequency analysis and processing. We selected the fourth-order daubechies wavelet as the mother wavelet that was similar to the waveform of the PA signal to perform 8-layer wavelet decomposition on the photoacoustic signals. Considering effectiveness and simplicity, the sqrtwolog threshold was chosen and the calculation formula is[38]

where M is the signal length and σ is the estimated value of noise variance. The threshold type is soft threshold and the threshold processing can be expressed as

where c and cnew are the wavelet coefficients before and after the threshold processing. In addition to white noise, we can also see clutters in Fig. 2(a). Wavelet analysis can effectively reduce white noise, but the signal is also disturbed by clutters. We selected the signal with better waveform which was not affected by clutters as the template and calculated the correlation coefficients between the collected signals and the template signal after wavelet denoising to remove the signals in which the clutters and the photoacoustic signal overlapped in time domain. The normalized correlation coefficient is given by[39]

where γxy is the normalized correlation coefficient, M is the length of the template signal, x(m) denotes the template point, y(m) denotes the signal point under analyzed, is the mean of the template points and is the mean of the signal points. Only the signal whose correlation coefficient was higher than the threshold value was retained. In our experiment, the correlation coefficient between the signal not disturbed by clutters and the template signal was usually above 0.8 while that between the signal disturbed by clutters and the template signal was usually below 0.4. The threshold value was set to 0.7, which could effectively distinguish the useful signals and the signals disturbed by clutters. In the processed signals like the red curve in Fig. 2(a), the white noise has been effectively reduced and the useful signals were not disturbed by clutters. The SNR of the denoised signals reached 11.06, which was 5.17 times higher than that before denoising. We effectively improved the SNR of the photoacoustic signals with wavelet denoising and correlation detection method while Wang et al. averaged 64 acquisitions of the photoacoustic signals but only raised the SNR to 3.9.[29] Sun et al. raised the SNR of the photoacoustic signals from 3.86 to 25.2.[31] In our experiment, the raw photoacoustic signals were weak and had low SNR which restricted the SNR of the processed signals.

As shown in Fig. 2(c), in the image reconstructed using the denoised signals, some sample points can effectively recognize while the others remain unrecognizable and the sample image is far from being recovering accurately. The raw photoacoustic signals were weak and prone to be affected by white noise and clutters. Although we collected 10 signals at each point, there were still many sample points where no useful signals were collected. This was the main reason why many points were still unrecognizable after processing. Obviously, when the intensity of photoacoustic signals is low, both signal extraction and scanning imaging are difficult.

In order to improve the photoacoustic signal amplitude effectively and obtain high quality photoacoustic image, we used DMD based superpixel method for wavefront shaping of the incident light. In wavefront shaping experiment, the enhancement[32] defined as the ratio of the intensity of the focus after and before optimization is often used to evaluate the effect of optimization. The theoretical enhancement of iterative optimization wavefront shaping can be expressed as

where N is the number of controllable degrees of freedom on the optical modulator and α is the relative enhancement decided by modulation mechanism. For pure amplitude modulation, α = 1/2π and for perfect amplitude and phase co-modulation, α = 1. DMD is a binary modulation digital micromirror device. In the superpixel method, DMD is divided into independent superpixel blocks. Each superpixel block is composed of n × n subpixels and each micromirror on DMD serves as a subpixel. The light reflected from the DMD surface passes through a 4f system consisting of two lenses and a spatial filter. The position of the spatial filter is specially selected that only the +1 or −1 diffraction beam could pass it through. The light reflected by different subpixels in the same superpixel block will arrive at the spatial filter with different phases. Phase of the superpixel block can be adjusted by changing the subpixel combination since the superpixel block is the superposition of the subpixels. In our experiment, a superpixel block was consisted of 4 × 4 subpixels and four phases 0, π/2, π, 3π/2 were optional. Larger the optical modes, easier to adjust the speckle grains to the comparable size of the absorption sample. In our experiment, each optical modes contained 20 × 20 micromirrors, that is, 5 × 5 identical superpixel blocks. And a superpixel mask was consisted of 13 × 24 optical modes. Further enlarging the optical modes would decrease the number of controllable degrees of freedom and decrease the enhancement. High robustness genetic algorithm[40] was used for iterative optimization and each iteration contained 20 masks. Twenty random superpixel masks were generated as the initial generation and the collected photoacoustic signals served as the feedback signal of the light intensity at the target point on the sample to sort the 20 masks. The higher the ranking, the higher the probability of being selected to generate the masks of the next generation through crossing and mutating. As the iterative optimization proceeded, the photoacoustic signals enhanced gradually and tended to convergence as shown in Fig. 3(a). In the experiment, we carried out 80 iterations. Further increase the number of iterations, the enhancement of the photoacoustic signal was not significant while the optimization time became longer. Figure 3(b) shows in detail the first generation, the 20th generation, the 40th generation, the 60th generation photoacoustic signals and the optimal photoacoustic signal.

Fig. 3. (a) Evolution curve of the averaged photoacoustic signal amplitude. (b) Photoacoustic signals collected during the wavefront shaping optimization. The black, red, blue, purple, and green curves represent photoacoustic signals of different generations in the optimization process respectively.

After optimization, the scattered light was effectively focused at the target point on the sample and the photoacoustic signal was greatly enhanced with the optimal mask loaded on the DMD as shown in Fig. 4. The black curve represents the random photoacoustic signal collected when a random superpixel mask was loaded on the DMD. The red curve represents the photoacoustic signal with amplitude 1 and phase 0 for every superpixel block which is referred to as the full white superpixel mask. The blue curve represents the optimal photoacoustic signal obtained after 80 iterations, which is 10.46 times higher than the random photoacoustic signal and 7.16 times higher than the photoacoustic signal of the full white superpixel mask. The transverse focal diameter of the focused UT can be calculated according to the theoretical formula[35]

where cs is the velocity of ultrasound in water, FUT is the focal length of the ultrasonic transducer, fUT is the central frequency of the ultrasonic transducer, and D is the crystal diameter of the ultrasonic transducer. The calculated transverse focal diameter of the 5-MHz UT is about 880 μm. The amount of the speckle grains on the sample in the detection region can be calculated as G = 880 μm × 150μm/[π × (75 μm)2] ≈ 7.47. The theoretical enhancement for pure amplitude modulation is η = 1/2π × 13× 24/7.47 ≈ 6.68 and that for perfect amplitude and phase co-modulation can be η = 1× 13× 24/7.47 ≈ 41.77. The experimental enhancement with DMD based superpixel method did not reach the theoretical value because we used only four phase steps in the process of amplitude and phase co-modulation of the incident light and the low signal-to-noise ratio of the initial signal had a certain impact on signal enhancement. Nevertheless, the enhancement in our experiment, 10.46, is much higher than the theoretical value of pure amplitude optimization. Wang et al. and Sun et al. performed pure amplitude optimization of the incident light with DMD[29,31] Wang et al. divided the DMD into 1024 individual blocks and the size of the speckle grains was similar to the detection region of the UT, i.e., G = 1. The theoretical enhancement of the photoacoustic signals was 163 and the experimental enhancement was about 14. Sun et al. divided the DMD into 18 × 32 blocks and G ≈ 7.47 in their experiment. They obtained an approximately 7.83 times enhancement of the photoacoustic signals while the theoretical enhancement was 12.3. Neither of the two experiments achieved the theoretical enhancement of pure amplitude modulation.

Fig. 4. Comparison of three photoacoustic signals. The black, red, and blue curves represent the photoacoustic signals generated when the random mask, the full white mask, and the optimal mask loaded on the DMD, respectively.

After optimization, the optimal mask was loaded on the DMD and remained unchanged. We moved the sample with the two-dimensional translation stage and collected the photoacoustic signals to reconstruct the image. Firstly, we used the collected photoacoustic signals to reconstruct the image directly without any processing on the signals. Figure 5(a) shows the scanning image reconstructed with the unprocessed signals. It is quite different from the scanning image before optimization. The absorption sample in the scanning image can be clearly distinguished and the image contrast reaches 0.84. Here, the contrast is defined as

named Weber contrast,[41] where I and Ib are the luminance of the sample image and the background. Figure 5(b) shows the image amplitude distribution along the twentieth column in Fig. 5(a) marked with the white line. The five peaks in Fig. 5(b) correspond to the five absorption nylon wires in Fig. 5(a). The half-width at those peaks are all no more than 3 points, corresponding to a resolution of at least 450 μm.

Fig. 5. (a) Scanning image with the unprocessed signals. (b) The image amplitude distribution along the white line in (a). (c) Scanning image with the denoised signals. (d) The image amplitude distribution along the white line in (c). The insets represent the optimal superpixel mask.

Then we processed the collected photoacoustic signals to remove the noise mixed in the signals. Figure 5(c) shows the scanning image reconstructed with the denoised signals. Compared with Fig. 5(a), the image contrast and resolution are significantly improved since the white noise and clutters mixed in the signals are removed. The image contrast reaches 8.57 due to the effective removal of background noise. Figure 5(d) shows the image amplitude distribution along the white line in Fig. 5(a). The half-width at the peaks in Fig. 5(d) are no more than 2 points, corresponding to a resolution of at least 300 μ m. Before optimization, the energy of the speckle field was dispersed and the photoacoustic signal was weak that the signal was almost submerged in white noise and clutters. It was difficult to extract the signals and image in this case. After optimization, the scattered light was focused at the target point on the sample and the photoacoustic signal significantly increased. Ultimately, we obtained much stronger photoacoustic signals and a high contrast photoacoustic image of the absorption sample through DMD based superpixel wavefront shaping method.

4. Conclusion

Photoacoustic imaging is challenged by weak photoacoustic signal and low signal-to-noise ratio. We proposed a new wavefront shaping and imaging method of low SNR photoacoustic signal based on DMD and superpixel method.

When the incident light was seriously affected by scattering, weak photoacoustic signals were excited and they were submerged in white noise and clutters. We performed wavelet analysis and correlation detection on the raw photoacoustic signals and raised their SNR from 2.14 to 11.06. We modulated the wavefront phase and amplitude of the incident light with DMD based superpixel method and the enhancement of the photoacoustic signal reached 10.46. The theoretical enhancement of pure amplitude modulation is only 6.68. Our experiment demonstrated that phase and amplitude co-modulation with DMD based superpixel method could achieve much better effect than pure amplitude modulation even when the SNR of photoacoustic signal was low. The enhancement can be further improved by raising the SNR of the original photoacoustic signals, applying more fine-grained phase control and utilizing the ultrasonic transducer with higher central frequency and tighter detection focal region.

After optimization, we kept the optimal mask loaded on the DMD and moved the absorption sample in the x–z plane by the two-dimensional translation stage to scan and image. A clear image with contrast of 8.57 and resolution of 300 μ m was obtained while imaging with original photoacoustic signals directly could not be achieved. The nylon wires could be clearly distinguished in the photoacoustic image. The imaging resolution was mainly limited by the ultrasonic transducer and better resolution can be achieved using ultrasonic transducer with higher central frequency.

Different from the time-consuming signal averaging method, the wavelet denoising and correlation detection of signals only took less than 50 ms. The optimization process took about 20 minutes and was mainly limited by the laser repetition rate of 10 Hz. At present, the refresh frequency of DMD can be as high as 32 kHz. With higher repetition rate pulse laser and faster calculating devices such as FPGA,[42] optimization can be accomplished in several seconds.

We experimentally demonstrated the superior performance of the DMD based superpixel method in photoacoustic wavefront shaping and imaging research of low SNR photoacoustic signals. Our research opens new perspectives for photoacoustic wavefront shaping and imaging with weak photoacoustic signals. The proposed method presents huge development potential and would promote the development of photoacoustic wavefront shaping technology and the practical application of photoacoustic wavefront shaping in biological imaging.

Reference
[1] Dravid A Mazumder N 2018 Laser Med. Sci. 33 1849
[2] Agarwal A Huang S W O’Donnell M Day K C Day M Kotov N Ashkenazi S 2007 J. Appl. Phys. 102 064701
[3] Majeed H Kandel M E Han K Luo Z Macias V Tangella K Balla A Popescu G 2015 J. Biomed. Opt. 20 111210
[4] Jacques S L 2013 Phys. Med. Biol. 58 R37
[5] Gareau D S Abeytunge S Rajadhyaksha M 2009 Opt. Lett. 34 3235
[6] Huang D Swanson E A Lin C P Schuman J S Stinson W G Chang W Hee M R Flotte T Gregory K Puliafito C A Fujimoto J G 1991 Science 254 1178
[7] Beaurepaire E Oheim M Mertz J 2001 Opt. Commun. 188 25
[8] Ntziachristos V 2010 Nat. Methods 7 603
[9] Duck F A 2013 Physical properties of tissues: a comprehensive reference book London Academic Press 55 122 10.1016/B978-0-12-222800-1.50007-3
[10] Xu M Wang L V 2006 Rev. Sci. Instrum. 77 041101
[11] Yuan Y Yang S 2012 Chin. Phys. 21 054211
[12] Yu H Park J Lee K R Yoon J Kim K D Lee S Park Y K 2015 Curr. Appl. Phys. 15 632
[13] Popoff S M Lerosey G Carminati R Fink M Boccara A C Gigan S 2010 Phys. Rev. Lett. 104 100601
[14] Zhang Z Zhang B Feng Q He H Ding Y 2018 Chin. Phys. 27 084201
[15] Yu L Xu X Zhang Z Feng Q Zhang B Ding Y Liu Q 2019 Chin. Phys. Lett. 36 114203
[16] Vellekoop I M Mosk A P 2008 Opt. Commun. 281 3071
[17] Cui M Yang C 2010 Opt. Express 18 3444
[18] Shen Y Liu Y Ma C Wang L V 2016 J. Biomed. Opt. 21 085001
[19] Xu X Liu H Wang L V 2011 Nat. Photon. 5 154
[20] Yu Z Li H Lai P 2017 Appl. Sci. 7 1320
[21] Kong F Silverman R H Liu L Chitnis P V Lee K K Chen Y C 2011 Opt. Lett. 36 2053
[22] Chaigne T Katz O Boccara A C Fink M Bossy E Gigan S 2014 Nat. Photon. 8 58
[23] Caravaca-Aguirre A M Conkey D B Dove J D Ju H Murray T W Piestun R 2013 Opt. Express 21 26671
[24] Conkey D B Caravaca-Aguirre A M Dove J D Ju H Murray T W Piestun R 2015 Nat. Commun. 6 7902
[25] Chaigne T Gateau J Katz O Bossy E Gigan S 2014 Opt. Lett. 39 2664
[26] Cao Z Mu Q Hu L Lui Y Xuan L 2007 Chin. Phys. 16 1665
[27] Jang M Ruan H Vellekoop I M Judkewitz B Chung E Yang C 2015 Biomed. Opt. Express 6 72
[28] Gu C Zhang D Chang Y Chen S C 2015 Opt. Lett. 40 2870
[29] Tay J W Liang J Wang L V 2014 Opt. Lett. 39 5499
[30] Tzang O Piestun R 2016 Opt. Express 24 28122
[31] Sun J Zhang B Feng Q He H Ding Y Liu Q 2019 Sci. Rep. 9 4328
[32] Vellekoop I M 2015 Opt. Express 23 12189
[33] Goorden S A Bertolotti J Mosk A P 2014 Opt. Express 22 17999
[34] Jia Y Feng Q Zhang B Wang W Lin C Ding Y 2018 Chin. Phys. Lett. 35 054203
[35] Chaigne T 2016 Control of scattered coherent light and photoacoustic imaging: toward light focusing in deep tissue and enhanced, subacoustic resolution photoacoustic imaging Ph.D. Dissertation Paris Pierre and Marie Curie University https://www.theses.fr/2016PA066162
[36] Qiu H Lee J Lin J Yu G 2006 J. Sound Vib. 289 1066
[37] Peng Z Wang G 2017 Sci. Rep. 7 4564
[38] Donoho D L Johnstone I M 1994 Biometrika 81 425
[39] Norollah M Pourghassem H Mahdavi-nasab H 2012 Fourth International Conference on Computational Intelligence and Communication Networks November 3–5, 2012 Mathura, India 331 10.1109/CICN.2012.129
[40] Conkey D B Brown A N Caravaca-Aguirre A M Piestun R 2012 Opt. Express 20 4840
[41] Li W Lu C Zhang J 2013 Opt. Laser Technol. 45 654
[42] Wang D Zhou E H Brake J Ruan H Jang M Yang C 2015 Optica 2 728